scratchpadfandomcom_ja-20200214-history
PMT: V: SrTiO3 tutorial
Thanks to the robustness in the MTO+APW(PMT) method, it is possible to set up a ctrl file easily. Here are minimum instructions to set up a ctrl file for MTO+APW. (Read also MarksOriginalDoc/fp.html. The floating orbitals or orbital optimization are unnecessary for MTO+APW.) To set up a ctrl file, you first have to write down crystal structure in the ctrls.* file, e.g. ctrls.srtio3. Then you can run ctrlgen.py, e.g., "ctrlgen.py srtio3". ctrlgen.py convert ctrls.srtio3 into ctrl.srtio3, which can be used by the following lmf calculaiton (do lmfa in advance). We explain these procedure below. We expect "ctrlgen.py" should work well, and we will improve it in future. But it is not well tested yet. So let me know something strange happens. =How to make a ctrl file = All the files in this section is in TESTsamples/SrTiO3_tutorial. Prepare crystal structure file. "ctrls.{name}" Let's make a file which contains only a crystal structure. %const da=0 au=0.529177 #<-- These % sections are to define variables. %const d0=1.95/au a0=2*d0 v0=a0^3 vfrac=1 v=v0*vfrac a1=v^(1/3) % show vars HEADER SrTiO3 cubic #<-- Any header STRUC ALAT={a1} # NBAS: number of site, NSPEC: number of atomic species DALAT={da} PLAT=1 0 0 0 1 0 0 0 1 # PLAT: primitive vector in ALAT. ALAT (in a.u.) SITE ATOM=Sr POS=1/2 1/2 1/2 # Atom SPEC and position in cartesian coordinate ATOM=Ti POS= 0 0 0+0 ATOM=O POS=1/2 0 0 ATOM=O POS= 0 1/2 0 ATOM=O POS= 0 0 1/2 SPEC ATOM=Sr Z=38 ATOM=Ti Z=22 ATOM=O Z=8 Save it as ctrls.srtio3 :NOTE. :1. "% const" section defines variables, "% show vars" requests programs(lmf, lmfa, lmchk) to show defined variables at the begining of console output. As this example shows, defined variables are used as {something}, e.g, {a1}. :By command line arguments when you invoke lmf, it is possible to replace values of variables defined in the % const section by "--v{variable}={value}" argument when you run lmf or so. For example, " lmchk --va0=2.5 srtio3" uses a0=2.5 instead of a0=2*d0 defined in the % const section. (When you use this, this is recorded in save.* file in the case of lmf. See Li_Lattice/ sample and examine it; the lattice constant is reecorded in save.* file.) :2.The lattice constant is given by ALAT+DALAT. The computational cutoffs internally used within lmf is determined by ALAT only. In order to obtain total energy as a smooth function of lattice constant, it is necessary to to use DALAT when you change lattice constant. Generally speaking, there are ways to deform cells with keeping such cutoffs for crystal structure without deformation.Do lmf --input :Remember The tag for category, e.g, SPEC must start at the beginning of the line. The tags for tokens, e.g, ATOM=, Z= must not start at the beginning of the line. :3. STRUC_NBAS= and STRUC_NSPEC= are not necessary. Furthermore, you don't need to supply SPEC if you use standard name for atoms defined in "standard SPEC", which is shown when you run ctrlgen.py without arguments. Run a script ctrlgen.py Then run ctrlgen.py srtio3 It calls other programs. If an error occurs, probably your PATH is wrong. This ctrlgen.py is a python srcipt. It generates ctrl.srtio3 file (ctrls.srtio3 is unchanged). Here we explain how ctrlgen.py works (you can skip this explaination below). ctrl.gen.py Generate MT radius Internally(L733:ctrlgen.py) it calls "lmchk --getwsr tmp >llmchk_getwsr". (ctrl.tmp is ctrls.srtio3 + some default section). This generates suggested MT radius which are stored in rmax.tmp file. See fp.html and lmto.html#section8 for discussion about how to determine R. ("lmchk srtio3" shows geometrical informations of its crystal structure). Default setting of RSMH, EH, KMXA, PZ, P The ctrlgen.py adds these obtained R to a ctrl file named as ctrl.tmp2. Then it does "lmfa srtio3 >llmfa.tmp2" (L782:ctrlgen.py). It generates 'mtopara.tmp2', which contains suggested values for RSMH, EH, KMXA, PZ and P as Sr@ RSMH= 2.411 2.411 2.094 EH= -0.100 -0.100 -0.100 KMXA=7 PZ=0,14.9 P=0,5.3 Ti@ RSMH= 1.393 1.393 1.055 EH= -0.100 -0.100 -0.100 KMXA=7 O@ RSMH= 0.900 0.900 EH= -1.360 -0.512 KMXA=7 This is used as a part of ctrl file. Note: Console output (llmfa.tmp2) includes "conf: sections"; it shows atomic configurations for each SPEC_ATOM. In the case of srtio3, for Sr, it is conf:SPEC_ATOM= Sr : --- Table for atomic configuration --- conf int(P)z = int(P) where P is replaced by PZ if it is semicore conf: isp l int(P) int(P)z Qval Qcore CoreConf conf: 1 0 5 5 2.000 8.000 => 1,2,3,4, conf: 1 1 5 5 0.000 18.000 => 2,3,4, conf: 1 2 4 4 0.000 10.000 => 3, conf: 1 3 4 4 0.000 0.000 => conf: 1 4 5 5 0.000 0.000 => conf:----------------------------------------------------- (you can see this by "grep conf llmfa.tmp2".) Here int(P) means the principle quantum number for valence orbitals. (P is the continuous principle quantum number; integer part corresponds to the usual principle quantum number. See MarksOriginalDoc/lmto.html. A pair of RSMH and EH specify an envelope function (the smooth Hankel function). See fp.html. These PZ and P means that local orbital is with P=4.9, and the valence orbital is with 5.3 for p channel (first zero means no local orbital for s channel). 1 of 10th ditit in PZ means the local orbital is the extended local orbital (Logically speaking, it is not a local orbital but a MTO because it is extended and augmented). This default setting of mtopara.srtio3 is given in a subroutine writebasis in subs/ioorpb.F. It uses a default setting setting: Possible semicore is (Takao thinks) ----- these are important --- Zn,Ga,Ge : 3d PZ=0,0,13.9 P=0,0,4.2 Cd,In,Sn,Te : 4d PZ=0,0,14.9 P=0,0,5.2 Hg,Tl,Pb,Bi : 5d PZ=0,0,15.9 P=0,0,6.2 Lu,Hf,Ta,W : 4f PZ=0,0,0,14.9 P=0,0,0,5.15 somehow meaningful, but can be negligible --- Na,Mg : 2p PZ=0,12.9 P=0,0,13.2 K,Ca : 3p PZ=0,13.9 P=0,0,14.2 Rb,Sr : 4p PZ=0,14.9 P=0,0,15.2 Cs,Ba : 5p PZ=0,15.9 P=0,0,16.2 --- less important --- Li : 1s PZ=11.9 P=2.6 Q=0,1 (this is commented as #PZ. In this case Q is needed. See explanation shown by lmfa without Q, and SPEC_ATOM_Q part shown by "lmfa --input") Finally, ctrlgen.py add some default categories to ctrl.srtio3. Look into it. We have to edit it by hand as follows. Edit the ctrl file, if neccesary This file contains all the parameters to run lmfa and lmf. The ctrl.srtio3 generated by ctrlgen.py from ctrls.srtio3 is given as # Do lmf --input to see all effective category and token # VERS LM=7 FP=7 # version check. Fixed. IO SHOW=T VERBOS=35 # SHOW=T shows readin data (and default setting at the begining of console output) # It is useful to check ctrl is read in correctly or not (equivalent with --show option). # # lerger VERBOSE gives more detailed console output. SYMGRP find # 'find' evaluate space-group symmetry automatically. # Usually 'find is OK', but lmf may use lower symmetry # because of numerical problem. # Do lmchk to check how it is evaluated. %const kmxa=7 # kmxa=7 is good for pwemax=5 or lower. # larger kmxa is better but time-consuming (maybe not the critical part for large systems). %const da=0 au=0.529177 %const d0=1.95/au a0=2*d0 v0=a0^3 vfrac=1 v=v0*vfrac a1=v^(1/3) % show vars HEADER SrTiO3 cubic STRUC ALAT={a1} DALAT={da} PLAT=1 0 0 0 1 0 0 0 1 NBAS= 5 NSPEC=3 SITE ATOM=Sr POS=1/2 1/2 1/2 ATOM=Ti POS= 0 0 0+0 ATOM=O POS=1/2 0 0 ATOM=O POS= 0 1/2 0 ATOM=O POS= 0 0 1/2 SPEC ATOM= Sr Z= 38 R=3.616323 RSMH= 2.411 2.411 2.094 EH= -0.100 -0.100 -0.100 PZ=0 14.9 P=0 5.3 KMXA={kmxa} LMXA=4 MMOM=0 0 0 0 ATOM= Ti Z= 22 R=2.089960 RSMH= 1.393 1.393 1.055 EH= -0.100 -0.100 -0.100 KMXA={kmxa} LMXA=4 MMOM=0 0 0 0 ATOM= O Z= 8 R=1.595007 RSMH= 0.900 0.900 EH= -1.360 -0.512 KMXA={kmxa} LMXA=4 MMOM=0 0 0 % const pwemax=3 nk=2 BZ NKABC={nk} {nk} {nk} # division of BZ for q points. METAL=3 # METAL=3 is safe setting. For insulator, METAL=0 is good enough. # When you plot dos, set SAVDOS=T and METAL=3, and with DOS=-1 1 (range) NPTS=2001 (division) even for insulator. # (SAVDOS, DOS, NPTS gives no side effect for self-consitency calculaiton). BZJOB=0 # BZJOB=0 (including Gamma point) or =1 (not including Gamma point). # In cases , BZJOB=1 makes calculation efficient. ITER CONV=1e-6 CONVC=1e-6 NIT=30 # An other choice is # ITER MIX=A2,b=.5,n=3 CONV=1e-6 CONVC=1e-6 NIT=20 # Practically results are independent from mixing procedure. HAM NSPIN=1 # Set NSPIN=2 for spin-polarize case; then set SPEC_MMOM (initial guess of magnetic polarization). FORCES=0 # 0: no force calculation, 1: forces calculaiton GMAX=9 # this is for real space mesh. See GetStarted. REL=T # T:Scaler relativistic, F:non rela. XCFUN=1 # =1 for Vosko-Ceperly-Alder; GGA is not yet. # XCFUN=2 shows a bug for Hydrogen atom. # (subs/evxc.F works only for XCFUN=1 if rho(up)=0 or rho(down)=0). PWMODE=11 # 10: MTO basis only (LMTO) PW basis is not used. # 11: APW+MTO (PMT) # 12: APW basis only (LAPW) MTO basis is not used. PWEMAX={pwemax} # (in Ry). When you use larger pwemax more than 5, be careful # about overcompleteness. See GetStarted. ELIND=-1 # this is only for accelaration of convergence. Not need to change. OPTIONS PFLOAT=1 # Q=band (this is quit switch if you like to add) NOTE 1.(again): The categories,e.g.,SPEC must start at the beginning of the line. The tokens,e.g, ATOM=, Z= must not start at the beginning of the line. NOTE 2. In each line, columns after "#" is neglected (these are comments). You may need to edit this ctrl file. There are other category-token to control lmf in detail (do lmf --input to see it). ctrlgen.py only inlcudes minimum things. NOTE: To check that lmf reads parameters correctly, do lmf --show srtio3. At the beginning of console output, supplied input is shown. * For metal, we need large NKABC=n1 n2 n3 * KMXA=7 is safe setting up to PWEMAX=5 (PW cutoff by 5Ry), as far as we tested. In cases you can use smaller number to reduce the computational time. * Mixing parametes ITER_MIX is a little difficult to understand. lmf chooses a mixing scheme automatically. But you can set it by yourself, e.g. as "ITER MIX=A2,b=.5,n=3". * For the magnetic case like Fe, you have to use HAM_NSPIN=2 (this expression means NSPIN=1 in HAM category). Further, add initial magnetic momemnt specification MMOM=0,0,2 (this means polarized in d channel). We now have complete ctrl file. You can run lmfa and lmf (repeat lmfa for safe; when PZ is added, need to generate density PZ as valence). * VERBOSE=35 or 40 might be better to observe all eigenvalues. * To start calculations, we recommend HAM_PWMODE=11 HAM_PWEMAX=3 After it is converged, you may enlarge PWEMAX to 4, 5 or so. Note that larger PWEMAX requires larger KMXA. For example, we have to use KMXA.ge.6 for SrTiO3 for PWEMAX=5. In this case, you will see only a little difference even when you use KMXA=20 (this is a convergence check about KMXA for PWEMAX=5). When we used KMXA <7, the calculation stoped suddenly during iteration cycle; it gives very unstable divergent behevior (out of convergence path; then the Hamiltonan becomes very strange; we will have to improve KMAX-related pat ---> "Gaussian Polynomial" expansion of smooth Hankel function at atomic site). *FTMESH defines real-space mesh for charge density. See a section below for setting. Instead of FTMESH, you can use GMAX to specify the real-space mesh. Usually GMAX=9 gives almost converged results. See HAM_FTMESH section below. Run lmf Now you are ready to run lmf. Do lmfa srtio3 Then you need to remove mix.* (for mixing data in previous itteration), moms.*(momentum data in previous itteration), rst.* (this contains self-consistent results) if you want to start over. If they exist, they are read by lmf. (This is a defalut setting. To change this beheviror, use --rs{something} options; do "lmf --help" and read MarksOriginalDoc/fp.html). lmf srtio3 OPTIONS_PFLOAT You have to use OPTIONS_PFLOAT=1 (bug fix of old versions). Without it, radial functions are not calculated at the center of gravity of the occupied states. (But the previous version somehow worked well). OPTIONS_PFLOAT=1 is usually stable and gives more rapid convergence than old versions. How to determine l-cutoff? (STRUC_NL; SPEC_ATOM_LMXA,LMX,LMXL) There are some kinds of l-cutoffs. As "lmf --input" shows, they are :SPEC_ATOM_LMXA: l-cutoff for augmentation. LMXA=NL+1 for defaults. :SPEC_ATOM_LMXL: lmax for which to accumulate rho,V in sphere :SPEC_ATOM_LMX : l-cutoff for basis (head part specfied. this corresponds to # of parameters in EH= (the same that for RSMH=). :STRUC_NL : global default lmax+1 for basis and augmentation Specification by SPEC_ATOM is stronger than the global default by STRUC_NL. This default is not effective if SPEC_ATOM_LMXA are spefcified for all the SPECs. Usually you neither need to set LMX nor LMXL because they are automatically choosen in defaults (LMX= is given by the numbers of EH=,RSMH= sections. LMXL is chosed to be LMXA). I(takao) think a possible (automatic) choice of LMXA is to set LMXA= {# of EH channel} + 2 (or +1 ?) in MTO+APW mode. This is given by ctrlgen.py. Need to be tested. Read MarkOriginalDoc/fp.html to know what they mean. How to determine real-space mesh (HAM_FTMESH or HAM_GMAX). The smooth part of electron density is represented on a real-space mesh. "HAM_FTMESH {3 integers}" is the number of divisions along the three lattice vectors. Or you can use HAM_GMAX, which is the maximum size of wave vector. (if GMAX exists, FTMESH is neglected). Search these words in MarksOriginalDoc/fp.html. FTMESH is useful to set number of real-space mesh explicitly (when you require smoothness of total energy when you change lattice constant). Too small FTMESH will not give converged total energy. If FTMESH is too large, however, current version causes a numerical problem. It will be fixed in future. Thus, you may need to observe the change of the total energy (and eigenvalues) as function of FTMESH. Look into examples. In anyway the convergence as a function of FTMESH is rapid. There is another experimental observation to choose FTMESH in a reasonable manner; when you run lmf (run and stop it within seconds), console output near its begining shows that sugcut: make orbital-dependent reciprocal vector cutoffs for tol= 1.00E-05 spec l rsm eh gmax last term cutoff La 0 2.00 -0.10 3.393 1.16E-05 1055 La 1 1.25 -1.14 5.825 1.00E-05 5343 La 2 1.68 -0.10 4.551 1.01E-05 2509 La 3 1.21 -0.10 6.862 1.00E-05 8617 Ga 0 1.90 -0.55 3.574 1.01E-05 1237 Ga 1 1.86 -0.10 3.854 1.00E-05 1511 Ga 2 0.73 -0.10 11.052 1.01E-05 36243 O 0 0.88 -1.32 7.747 1.00E-05 12461 O 1 0.83 -0.38 8.930 1.00E-05 19109 look into "last term", you may need to tune GMAX or FTMESH so that "last term" is ~1e-05. (When "last term"is ~1d-5 or below, it indicates that the real-space mesh is dense enough, that calculations are converged well as for the mesh.) Overcompleteness When you use large PWEMAX (more than ~5Ry or so), be careful! It can cause a problem of poor linear-dependency of the basis sets; the overlap matrix of a basis set can have very small eigenvalues. The eigenfunction of the overlap matrix corresponding to such a small eigenvalue can be very strange; it is the linear combination of MTOs and APWs where there exist strong cancellations each other. confuse eigenvalues of the overlap matrix with those of the secular equation (eigenvalues of the Kohn-Sham equations). In such a case, it is necessary to reduce the dimension of the Hamiltonian (the Hilbert space spanned by MTOs+APWs) by eliminating some subspace (or removing some of basis functions). Without such a procedure, calculation stops suddenly, or the final results show very strange eigenvalues and total energies. There are two way (A) and (B) as explained below to overcome this linear dependency problem (OVNCUT is another way, but not explained here). HAM_OVEPS A way is to set HAM_OVEPS=1d-6 (or so). If # is a positive number, lmf will diagonalize the overlap matrix, and eliminate the subspace which has eigenvalues smaller than #. Note: In lm-7.0betaK, the space of MTO is preserved, and within the space spanned by APWs orthogonalized to the space of MTO, we calculate eigenvalues of overlap matrix. (See subroutine zhev_tk in slatsm/zhev.F). This procedure is a little different from lm-7.0beta (it does not preserve space of MTO). See fp.html. In the console output of lmf, "zhev_tk: ovlmat=" shows the eigenvalues of the overlap matrix. E.g., Li/llmf, which is an output of "lmf li", shows zhev_tk: ovlmat= 1 0.66D-05 2 0.14D-01 3 0.14D-01 4 0.14D-01 5 0.33D+00 6 0.10D+01 7 0.10D+01 8 0.10D+01 9 0.10D+01 10 0.10D+01 11 0.10D+01 12 0.10D+01 13 0.10D+01 14 0.92D+04 15 0.10D+05 16 0.10D+05 17 0.10D+05 18 0.11D+05 If we set HAM_OVEPS=1d-6, an eigenfunction corresponding to the eigenvalue 0.66D-05 of the overlap matrix is removed from the Hilbert space. Some of eigenvalues are ~1, this means these corresponds to APW which are linearly-independent. Other basis with 1D+5 are MTOs (when we calculate overlap matrix in zhev_tk, we set normalization (diagonal part of overlap matrix) as =1, and =1d+5, where {APW(q+G)} denotes the set of the APW, and {MTO_i(q)} denotes that of the MTO. With this weighted procedure, we can easily preserve the Hilbert space of MTO). Remove MTOs for the channel Removing some of MTOs explicitly can also reduce the Hilbert space. This method is better than (A) in practice since (A) can cause a sudden change of total energy when you change lattice constant(s). To remove the "Li 2s", set "RSMH=0 1.6 EH=0 -0.1" instead of original setting as "RSMH=1.6 1.6 EH=-0.1 -0.1". Then the smallest eigenvalue 0.66D-05 disappears. This is because the very small eigenvalue 0.66D-05 originates from the fact that the MTO of Li 2s is already well represented by the superposition of APWs. :remove MTOs of a channel, set zero for RSMH and EH for the channel as we explained above. Note: To remove all MTO from an atom use many zero (number of LMXA) as RSMH= 0 0 0 0 0 EH= 0 0 0 0 0 ; this is a minor bug. Checking convergence It is better to watch all the Kohn-Sham eigenvalues (shown after "eigenvalue=" in console output;) at least at the Gamma point or so; some of eigenvalues from the top becomes too high when the linear-dependency becomes poorer; however, the convergence of lmf can be not so bad as far as it converges since the self-consistency in LDA is only related to the eigenfunctions below the Fermi energy. Generally speaking, PWEMAX=5 with minimum MTOs gives not so bad result as for total energy (the energy which can be used to determine lattice constant or so; but the numerical convergence for each systems are dependent on sysmtems). When you use PWEMAX=10 or more, you must need to use procedure (A) or (B). Don't forget to use large KMXA (KMXA=15 or so is needed for PWEMAX=10) (note: to start over, remove mix* moms* rst*, and repeat it). How to relax atomic positions For relaxation of atoms (atomic force is calculated and relaxed), see a sample, LaGaO_323. This is a perovskite with twenty atoms per cell. What you have to do is to add HAM FORCES=1 DYN MSTATHESS=t XTOL=.001 GTOL=0 STEP=.015 NIT=30 "grep ATOM log.lagao" ---> it shows how the atomic positions are moved. You may need to add HAM_FRZWF=1 to calculate better force.